#!/bin/bash -e

function info() {
echo Usage: `basename $0` 'bam'
exit 1
}

while getopts  ":p:f:i:" opt; do
	case  $opt  in
		p) out_prefix=$OPTARG;;
		f) suffix=$OPTARG;;
        i) interval=$OPTARG;;
		*) info;;
	esac
done
shift $(($OPTIND - 1))


if [ $# -lt 1 ]; then info; fi

. $var

if test -n "$interval"; then interval0="-l $interval"; fi
if test -n "$interval"; then interval="-L $interval"; fi


$java_run/gatk38 \
	-T UnifiedGenotyper \
	-R $ref_genome $interval \
	-I $1 \
	-o $out_prefix.snp.vcf \
	-glm SNP \
	--downsample_to_coverage 20000 \
	-nt 8 \
	-stand_call_conf 0 \
	--min_base_quality_score 0 \
	--max_deletion_fraction 2 \
	-slod \
	--useNewAFCalculator \
	-hets 0.001 \
	--output_mode EMIT_ALL_SITES \
	--allSitePLs
	# --disable_read_filter

~/../fengbo.zeng/local/bin/samtools mpileup -u -v -f $ref_genome -d 20000 $interval0 -t DP,AD,ADF,ADR,SP $1 > $out_prefix.mpileup.vcf
bcftools call -Ov --multiallelic-caller --keep-alts -fGQ,GP $out_prefix.mpileup.vcf > $out_prefix.bcf.vcf

:<<mark
$java_run/gatk38 \
	-T UnifiedGenotyper \
	-R $ref_genome $interval \
	-I $1 \
	-o $out_prefix.indel.vcf \
	-glm INDEL \
	--baqGapOpenPenalty 45.0 \
	--downsample_to_coverage 20000 \
	-nt 8 \
	--min_indel_fraction_per_sample 0.001 \
	--min_indel_count_for_genotyping 1 \
	-stand_call_conf 0.1 \
	--min_base_quality_score 17 \
	--max_deletion_fraction 2 \
	-slod \
	--useNewAFCalculator \
	--output_mode EMIT_ALL_SITES

$java_run/gatk38 \
	-T UnifiedGenotyper \
	-R $ref_genome $interval \
	-I $1 \
	-o $out_prefix.indel.pairhmm.vcf \
	-glm INDEL \
	--baqGapOpenPenalty 40.0 \
	--downsample_to_coverage 20000 \
	-nt 8 \
	--min_indel_fraction_per_sample 0.001 \
	--min_indel_count_for_genotyping 1 \
	-stand_call_conf 0.1 \
	--min_base_quality_score 17 \
	--max_deletion_fraction 2 \
	-pairHMM LOGLESS_CACHING\
	-slod \
	--useNewAFCalculator

mark
. $cmd_done
